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ABSTRACT 

Simulations of core-collapse supernovae (CCSNe) result in successful explosions once the neutrino lumi- 
nosity exceeds a critical curve, and recent simulations indicate that turbulence further enables explosion by 
reducing this critical neutrino luminosity. We propose a theoretical framework to derive this result and take the 
first steps by deriving the governing mean-field equations. Using Reynolds decomposition, we decompose flow 
variables into background and turbulent flows and derive self-consistent averaged equations for their evolution. 
As basic requirements for the CCSN problem, these equations naturally incorporate steady-state accretion, 
neutrino heating and cooling, non-zero entropy gradients, and turbulence terms associated with buoyant driv- 
ing, redistribution, and dissipation. Furthermore, analysis of two-dimensional (2D) CCSN simulations validate 
these Reynolds-averaged equations, and we show that the physics of turbulence entirely accounts for the differ- 
ences between ID and 2D CCSN simulations. As a prelude to deriving the reduction in the critical luminosity, 
we identify the turbulent terms that most influence the conditions for explosion. Generically, turbulence equa- 
tions require closure models, but these closure models depend upon the macroscopic properties of the flow. To 
derive a closure model that is appropriate for CCSNe, we cull the literature for relevant closure models and 
compare each with 2D simulations. These models employ local closure approximations and fail to reproduce 
the global properties of neutrino-driven turbulence. Motivated by the generic failure of these local models, we 
propose an original model for turbulence which incorporates global properties of the flow. This global model 
accurately reproduces the turbulence profiles and evolution of 2D CCSN simulations. 

Subject headings: convection — hydrodynamics — instabilities — methods: analytical — methods: numerical 
— shock waves — supernovae: general — turbulence 



1. INTRODUCTION 

Though our understanding of the core-collapse super- 
nova (CCSN) explosion mechanism remains incomplete, re- 
cent simulations indicate that it is likely to involve multi- 
dimensional effects. In fact, in all proposed mechanisms 
neutrino-driven convection plays an important, if not vital, 
role. Motivated by these results, we present a theoretical 
framework to investigate the role of turbulence in launching 
successful explosions. Furthermore, we lay down the foun- 
dation for this framework by deriving self-consistent steady- 
state equations for the background and turbulent flows. 

The fundamental problem of CCSN theory is to deter- 
mine how the stalled shock transitions into a dynamic ex- 
plosion. Within a few milliseconds after bounce, the core- 
bounce shock wave stalls into an accretion shock dMazurekl 
Il982t iBruennll 19851 119891) . Unchecked, co ntinued accretion 
throu gh the shock would form a black hole dO' Connor & Ottl 
1201 ll). However, the preponderance of observed neut ron stars 
(iLorimeT et al.l2006l) and supernova (SN) explosions dLi et al.l 
dictates that the stalled shock is revived into an explo- 
sion most of the time. 

For more than two decades, the favored mechanism for core 
colla pse supernovae has be en the delayed-neutrino mecha- 
nism dBethe & Wilsonl[l985l) . In this model, a neutrino lumi- 
nosity of several times 10 52 erg/s cools the protoneutron star 
and heats the region below the shock. Under the correct con- 
ditions, this heating by neutrinos can revive the shock and pro- 
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duce an explosion. Unfortunately, except for the least massive 
stars which undergo core collapse, most of the detailed ID 
neutrino-transport-hydrodynami c simulations do not produce 
solutions containing ex p losions (ILiebendorfer et al.ll2.001blah 

Thompson et al. 



Buras et al.l 120031, 



120051: iKitaura et alil2006l) . 



Rampp & Jankal 120021: 
20031: ILiebendorfer et al 

While most ID simulations fail to produce explo- 
sive solutions, recent 2D simulations show promising 
trends. These simulations capture multi-dimensional in- 
stabilities th at aide the neu tr ino m e chanism in driving 
explos ion s dMarek & Jankal 120071 iMurphv & Burrows! 



2008bt iScheck et al.l 1 20081 iMarek & Jankal 



2009 



Fernandez & Thompson! 120091: ISuwa et al.l 12010 : 
Nordhaus et al 1 120101) . These instabilities include neutrino- 
driven convection (IBurrows et all 119951: iJanka & Miillerl 
119951: IMurphv & Burrows! I2008b[ iNordhaus et al.l 120101) 
and the standing a ccretion shock ins ta bility (SASI) 
dBlondin et al.l 120031: fOhnishi et .al.l 120061: IBurrows et all 
20061: IMarek & Jankal I2007L 120091: iFoglizzo & Tagger! 120001: 
Foglizzo 2009). Besides these two instabilities, other multi- 
dimensional processes may revive the s talled shock, including 
magnetohydrodyna mic (MHD) jets (IBurrows et al.l l2007at 
iDessart et al.l 120081) and acoustic power derived from asym- 
metric accretion and an oscillating PNS (IBurrows et al.ll2006l 
l2007bl) . Though the prevalence of these last two processes 
is sti ll debated (IDessart et al.l 120081: Weinberg & Ouataerj 
120081) . two points remain clear: (1) the solution to the CCSN 
problem is likely to depend on multi-dimensional effects, and 
(2) in all proposed mechanisms, neutrinos and turbulence 
play an important, if not central, role. Hence, whatever the 
mechanism, it is important to understand the role of neutrinos 
a nd turbulence. 
IBurrows & Goshyl (119931) proposed that if a critical neu- 
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trino luminosity is exceeded for a given mass accretion rate, 
the neutrino mechanism succeeds. These authors developed 
a steady-state model for an accretion shock in the presence 
of a parametrized neutrino heating and cooling profile. The 
two most important parameters of the steady-state model are 
the neutrino luminosity L Ve and mass accretion rate M. They 
found no steady-state solutions for luminosities above a criti- 
cal curve L v<!}Cr it = f(M), and interpreted this curve as separat- 
ing steady-state (or failed supernova) solutions from explosive 
solutions. However, this work did not prove that the solutions 
above the critical curve are in fact explosive, nor did they con- 
sider multidimensional effects. 

Using a similar neu trino parametrization as in the 
IBurrows & Goshvl d 19931) work, IMurphv & Burrows! d2008bl) 
showed in ID and 2D simulations that the solutions above 
the critical luminosity are in fact explosive. Moreover, they 
found that the critical luminosity in the 2D simulations is 
^70% of that in ID for a give n M. Additional investi- 
gations by Nordhaus et af| (120101) show that the critical lu- 
minosity is even further reduced in 3D. These results sug- 
gest that the critical luminosity is a useful theoretical frame- 
work for describing the conditions for successful ex plosions 
(IBurrows & Goshvil 19931: [Murphy & Burrowsll200 8b). 

Initial investigations by IMurphv & Burrows! d2008bl) also 
suggest that the reduction in the critical luminosity is caused 
by turbulence. An alternative but related condition for 
explosion i s a comparison of the advec t ion and heating 
timescales (iThompsonl 120001: iJankal 120011 : iThompson et alj 
12001 IMurphv & Burrowsll2008bl) . The heating timescale is 
the time it takes to significantly heat a parcel of matter, and 
the advection timescale is the time to advect through this 
region. If the advection timescale is long compared to the 
heating timescale, then explosion ensues. In ID, as mat- 
ter accretes onto the PNS, it is limited to advect through 
the gain region with one short timescale. In 2D, convec- 
tive motions increase the dwell time, which leads to more 
heating for the same neutrino luminos it y and a lower critical 
lumino sity dMurphy & Burrowsll2008bT) . iPejcha & Thompsonl 
d201 ll) have recently challenged this explanation and suggest 
that rather than increasing the heating, turbulence acts to re- 
duce the cooling. Regardless, the simulations show that the 
critical luminosity is lower in the presence of convection. 

These results suggest that a theory for successful explo- 
sions requires a theoretical framework for turbulence and its 
influence on the critical luminosity. In this paper, we develop 
the foundation for such a framework. Recent developments 
in turbulence theory have led to accurate turbulence models, 
and in this paper we use similar strategies to develop a turbu- 
lence model appropriate for CCSNe. Such a turbulence model 
can then be incorporated into steady-state accretion models 
to derive reduced critical luminosities for explosion, as well 
as used in ID radiation-hydrodynamic (RHD) simulations to 
expedite systematic studies of core-collapse physics. In the 
present paper, we develop a turbulence model which captures 
the salient features of 2D core collapse convection, but even- 
tually the model must be calibrated against 3D simulations. 

To develop a turbulence model appropriate for core- 
collapse turbulence, we use a general and fu lly self-consistent 
appro ach c alled Reynolds decomposit ion dPlate et al] 1 19971 : 
lPopdl2000l:lLaunder & Sandhamll2002l) . The first step in this 
approach is to decompose the flow variables into averaged and 
fluctuating components. Evolution equations for the mean 
flow variables are then developed by writing the conserva- 



tion laws in terms of these mean and fluctuating components 
and then averaging. The resulting evolution equations for the 
mean fields contain terms which involve both the mean fields 
as well as correlations of fluctuating components. These cor- 
relations represent the action of turbulence and include the 
Reynolds stress, turbulent kinetic energy, turbulent enthalpy 
flux, and higher order correlations. These evolution equa- 
tions are self-consistent and naturally include the effects of a 
background flow, which is important in core collapse. Unfor- 
tunately, this procedure always produces evolution equations 
which depend on a correlation of higher order than the evolu- 
tion equations. Therefore, in order to develop a closed system 
of equations, the highest order correlations must be modeled 
in terms of the lower order correlations and mean fields. Fur- 
thermore, closure depends upon the macroscopic flow itself, 
so there is no unique closure for turbulence. This is the infa- 
mous closure problem of turbulence. 

Fortunately, there is a small class of turbulent flows 
(e.g., shear, buoyancy driven, etc.) and closure models 
have been developed that work well for each type un- 
der a rang e of conditions dTurnerl 119731: iPlate et al.l 119971 : 
|R^|200(I ILaunder & Sandhaml 120021: 1 Wilcox! 120061) . The 
general strategy to find closure relations involves an in- 
terplay b etween theory, observations , and numerical sim- 
ulations dLaunder & Sandhaml 120021) . First, terms in 
the mean-field equations are compared to either obser- 
vations or numerical simulations. Approximations are 
then proposed for the higher order correlation terms that 
satisfy the observations/simulations and provide closure. 
This appr oach has been used succ essfully for geophysi- 
cal flows dLaunder & Sandhaml 120021) a nd is now being ap- 
plied to stellar structu re calculations dGaraud et alj 120101 : 
iMeakin & Arnettll2007l) . 

Following on these successes, we use this strategy to de- 
velop a turbulence model for the core-collapse problem in 
which buoyancy and a background accretion flow dominate. 
In §|2j we use Reynolds decomposition to formally derive the 
averaged background and turbulence equations and identify 
terms that are important for neutrino-driven convection. Us- 
ing 2D simulations, we examine in §[3]the turbulent properties 
of neutrino-driven convection and show that the turbulence 
equations which we derive in §|2]are consistent with the sim- 
ulated flows. Finding solutions to the mean-field equations 
requires a closure model. Therefore, in §|4]we present several 
models representative of the literature. However, these fail to 
reproduce the global profiles of neutrino-driven convection, 
leading us to develop a novel global model. In we com- 
pare the results of the turbulence models (§ HI with the re- 
sults of 2D simulations and conclude that our global model is 
the only model to reproduce the global properties of neutrino- 
driven convection. In §[6j using the mean-field equations and 
2D simulations, we investigate the effects of turbulence on 
the conditions for successful explosions. Finally, in $7] we 
summarize our findings, and motivate the need for a similar 
analysis using 3D simulation data. 

2. BACKGROUND AND TURBULENCE EQUATIONS: REYNOLDS 
DECOMPOSITION 

The first step in understanding the effects of turbulence is 
to derive the governing steady-state equations. Therefore, in 
this section, we use Reynolds decomposition to derive exact 
equations for the steady-state background and turbulent flows. 

Historically, turbulence modelers have used two approaches 
to derive the models. In one, ad hoc equations are sug- 
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gested to model the important turbulence physics. These 
are often of practical use, but the underlying assumptions 
often make these of limited use. Mi xing length theory 
MLT) is one such approach (see {34.51 for its limitations) 
Kippenhahn & Wiegerjll990l) . In the second approach, one 
derives self-consistent equations for turbulence by decompos- 
ing the hydrodynamic equations into background and turbu- 
lent parts. Reynolds decomposition is an example of this 
approach. Though these equations are exact, they are not 
complete; they need a model for closure. If the ad hoc ap- 
proaches accurately represent Nature, then one should be able 
to derive them by making the appropriate assumptions in the 
exact equations. Therefore, regardless of the technique em- 
ployed, starting with the self-consistent equations enables a 
better understanding of the assumptions and limitations. In 
this paper, we pursue both approaches, but in this section, we 
use Reynolds decomposition to derive and explore the self- 
consistent equations for the background and turbulent flows. 

In Reynolds decomposition, the hydrodynamic equ ations 
are d ecomposed into background and turbulent flows (Pope 
120001) . Consider a generic flow variable, </>, and its decompo- 
sition into average (background) and fluctuating (turbulent) 
components: <f) = (f>o + <fi' . The mean-field background of <fi, 
(4>), is obtained by coarse spatial and temporal averages. The 
interval for the averages must be large or long enough to 
smooth out short term turbulent fluctuations, but they must 
not be too large or long so that interesting spatial or tempo- 
ral trends in the mean-field quantities are completely averaged 
out. Choosing the scales of the averaging window is depen- 
dent upon the problem, and in this paper, we define () as aver- 
aging over the solid angle in the spherical coordinate system 
and over a fraction of the eddy crossing time of the convective 
region. By definition, the coarse average of (f> is ((f)) = <fro and 
the mean-field average of the fluctuation is identically zero, 
(4>') = 0. Therefore, first order moments of turbulent fluctua- 
tions are identically zero and only higher order terms survive. 
For example, the average of the velocity fluctuation is zero, 
(v'), but the mean-field of the second order term, the Reynolds 
stress (R) = (v' ® v'), is nonzero. 

2.1. Averaged Background Equations 

The general equations for mass, momentum, and entropy 
conservation are 

^ + V-(p«) = 0, (1) 

^^ + V-(/9M(8)M)=-VP + pg, (2) 

ot 

and 

r(^ + v.(p*))=tf. (3) 

In these equations, the density, velocity, pressure, tempera- 
ture, and specific entropy are p, u, P, T, and s. g is the gravi- 
tational acceleration, and q is the specific local heating and/or 
cooling rate. 

After Reynolds decomposition, averaging, and assuming 
steady state, the hydrodynamics equations, eqs. ([HO, become 

V-( Po v+(pV)) = 0, (4) 

(pH)-Vv = -VPo + A#-V-(pR) (5) 

and 

<p«).V,o = (f) + ^-V.<F s >. (6) 



For all quantities, the turbulent perturbations are denoted by 
a superscript /, and, with the exception of velocity, the back- 
ground flow is denoted by subscript 0. The background ve- 
locity is v, and the perturbed velocity is v'. 

Equations @]|6] are very similar to the usual steady-state 
equations of hydrodynamics, but the last term in all three 
equations add new turbulence physics. Conservation of mass 
flux is split between the background and the turbulence, 
(p'v'). In the momentum equation, the extra force due to tur- 
bulence is the divergence of the Reynolds stress, R = v' ® v'. 
The entropy equation has two new terms; the divergence of 
the entropy flux, F s = ps'v', represents entropy redistribution 
by turbulence, and e represents heat due to turbulence dissipa- 
tion. For isotropic turbulence, the divergence of R can be re- 
cast as the gradient of turbulent pressure: i.e. -(2/3)V(pK), 
where K = v' ■ v'/2 is the turbulent kinetic energy. Using 
thermodynamic relations, we reduce the number of turbulent 
correlations in eqs. (@]|6]) by noting that (p'v 1 ) = r)F s , where 
r\ = pj/cp, cp is the specific heat at constant pressure, and 
$t = —(d]np/d]nT)p is the logarithmic derivative of density 
with respect to temperature at constant pressure. 

While we consider the convective entropy flux, F s , tra- 
ditionally, astrophysicists have considered the enthalpy flux, 
F e = pocp (v'T'} in turbulence models. The enthalpy flux has 
units of energy flux. Therefore, the enthalpy flux is a natural 
choice for stellar structure calculations in which the enthalpy 
flux and radiative flux must add to give the total luminosity 
of the star. Because the convective region is semi-transparent 
to neutrinos, there is no such constraint in the core-collapse 
problem. Furthermore, since we decompose the entropy equa- 
tion, the entropy flux is the most natural flux to consider. For 
the instances that require discussing the enthalpy flux, we ex- 
press it in terms of the entropy flux, F e = TqF s . 

In expressing {(p'v'} = r}F s , we have eliminated one of the 
turbulent correlations , but there still remain three turbulent 
correlations (R, F s , and e), resulting in more unknowns than 
equations. To clos e these equations, we derive the turbulence 
equations in § 12.2b . 

2.2. Averaged Equations for Turbulent Correlations 

Using the definitions of R and F s and the conserva- 
tion equations, we re-derive the evolution equations for the 
Reynolds stress (R) and the entropy fl ux (F s ) . For sim- 
ilar derivations of t hese equations, see ICanutol d 19931) and 
iGaraud et alj ( 120101) . The convective Reynolds stress equa- 
tion is 

d(pR)/dt + v-V(pR) + (pR)®V-v = 

+ (p'v') <E)g+ [(p'v') ®g] T Buoyant production 
- (pR) • Vv - [ (pR) ■ Vv] T Shear production 
-(V®F P +[X/®F P ] T ) Pressure flux (7) 

+(P'Vv'+[P'Vv'] T ) Pressure strain 

-V • (pv' <8 R) Turbulent transport 

-pos Dissipation , 

where the turbulent pressure flux is Fp = P'v', the dissipa- 
tion tensor is p e = p[V 2 (R) -2(Vv') • (Vv')], and [] T is the 
transpose operator. In eq. (0), we separate the terms on the 
right-hand-side into rows to better illustrate their physical rel- 
evance. They are buoyant and shear production, redistribution 
by the turbulent pressure flux, the pressure-strain correlation, 
turbulent Reynolds stress transport, and the turbulent dissipa- 
tion. In neutrino-driven convection of core collapse, buoyancy 
is the most important turbulent production. In terms of driv- 
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ing turbulence, the shear production term is less important. 
However, this term and the pressure-strain term are primarily 
responsible for redistributing stress among the components. 
For example, gravity acts mostly on the vertical stress com- 
ponents, but the shear production and pressure-strain terms 
redistribute stress to the horizontal components. Also impor- 
tant in redistributing stress is the turbulent transport term. In 
fact, in the next paragraph, we show that this term is in effect 
the divergence of the turbulent kinetic energy flux, which is 
very important in vertical kinetic energy transport. 

Taking the trace of the Reynolds stress equation gives the 
convective kinetic energy equation: 

d(pK) jdt + v ■ V (pK) + (pK) V • v = 

+ (pV)-g-tr((pR)-W) (8) 
-V-(F K }-V-(F P } + (P'V-v')-p e, 

where tr() is the trace operator, the turbulent dissipation be- 
comes e = tr(e)/2, and Fk = pv'K is the turbulent kinetic en- 
ergy flux. Once again on the right-hand-side, we have the fa- 
miliar terms: buoyancy and shear production, turbulent redis- 
tribution by the turbulent kinetic energy flux, the divergence 
of the pressure flux, work done by turbulent pressure, and tur- 
bulent dissipation. 
The corresponding equation for convective entropy flux is 

d{F s )/dt + v-V{F s ) + {F s )V-v = 

+PQT](Q)g Buoyant production 

- (F s ) ■ Vv- (pR-) Vio Gradient production 

- (s'VP') Pressure covariance ^ ' 
-V • (v' ®F S ) Turbulent transport 

+ (pv'q/T) heat production . 

where Q = s' 2 is the variance of the entropy perturbation. 
Once again, we separate the terms in eq. (0 into rows to high- 
light their physical significance. The terms in the entropy flux 
equation are analogous to the terms in the Reynolds stress 
equation. They are buoyant and gradient production terms, 
the pressure-covariance, turbulent transport, and heat produc- 
tion. Unlike the Reynolds stress equation, we find that buoy- 
ant production, gradient production, pressure covariance, and 
turbulent transport are all equally relevant in determining the 
entropy flux. 

The first term on the right hand side of eq. (0 is the buoy- 
ant production. This term is an important source in eq. ©, 
but it depends upon yet another correlation, the variance of 
the entropy perturbation. The corresponding equation for the 
entropy variance is 

d(pQ)/dt + V-(pv(Q)) = 

-2 (F s ) ■ Vso Gradient Production . . 

+2{s'Q/T) Heat Production 

-V • (pv'Q) Turbulent Transport . 

Equations dTKTOli are an exact set of evolution equations 
for the 2 nd order correlations (i.e. Reynolds stress, entropy 
flux, and entropy variance). While these equations are exact, 
they are not complete. Each equation depends upon 3 rd order 
correlations, necessitating further evolution equations for the 
higher order correlations. However, it is impossible to close 
the turbulence equations in this way, as each set of evolution 
equations depends upon yet higher order correlations. The 
only solution is to develop a closure model to relate higher 
order moments to lower order moments. 



This is analogous to the closure problem in deriving the 
hydrodynamics equations. In the hydrodynamics equations, 
the equation of state (EOS) is a microphysical closure model 
which relates the pressure (a higher moment) to the den- 
sity and internal energy. Because of the vast separation 
of scale, the EOS depends upon microphysical processes 
only and is independent of the macroscopic hydrodynam- 
ical flows. Hence, as a closure model, the EOS enables 
the hydrodynamic equations to be relevant for a wide range 
of macroscopic flows. Conversely, turbulence occupies the 
full range of scales from the microscopic to the largest bulk 
flows. In some cases, turbulence is the dominant macroscopic 
flow. Consequently, closure is necessarily dependent upon the 
macroscopic flow, making it impossible to derive a generic 
closure relation for turbulence. 

To find solutions to the turbulence equations, we need to 
construct a turbulence closure model that is appropriate for 
core collapse. The standard approach is to develop a turbu- 
lence closure model for each macroscopic flow. Fortunately, 
this task is not as daunting as it first appears. Turbulence can 
be divided into several classes that are characterized by the 
driving mechanism (i.e. shear, buoyancy, magnetic) and clo- 
sure models have been constructed that are appropriate for 
each class. For core collapse, buoyancy is the primary driving 
force, and the rest of this paper is devoted to finding an appro- 
priate buoyancy closure model for core-collapse turbulence. 

2.3. Steady-state Reynolds-averaged Equations in Spherical 

Symmetry 

Assuming a spherically symmetric background, the equa- 
tions for the background flow (eqs. [4]|6]l become 

M = 4tt r 2 m = 4tt r 2 (p v r + (p'v' r )) , (11) 

md r v r = -d r P Q + pogr - o r (pR rr ) + ^— , 

r r 

(12) 

and 

mdrSo = /^\ + ^-V r .(F r }, (13) 

where the r, 9, and <fr subscripts refer to the radial and angu- 
lar components in spherical coordinates. Therefore, d r is the 
partial derivative with respect to r, and V, • is the radial part 
of the divergence. Since we assume steady state, the mass 
accretion rate, M, is a constant. For isotropic turbulence, the 
last three terms of eq. (fT2l red uce to d r (\pK\ th e gradient 
of turbulent pressure. However. lArnett et alJ (120091) note that 
buoyancy-driven turbulence in spherical stars is not isotropic 
but is most consistent with R n = Ree+R<j><j> and Rgg = R<p<p. In 
this work, we adopt this later assumption where convenient, 
but we retain the general expression in eq. (fT2l) as a reminder 
that the relationships among the Reynolds stress components 
must be determined by theory, simulation, or experiment. 

The equivalent steady-state and spherically symmetric 
equations for Reynolds stress, entropy flux, and entropy vari- 
ation are 

V r d r (pR rr ) + (pRrr) V r ■ V r = 

+2(p'v' r )gr-2(pRrr)d r V r (14) 
- {v' r d r P' + d r P'V T ) - V r ■ (pv' r R rr ) - p e rr , 

v r d r {pR 9S )+ (pRee) V,- ■ v r = -V r • (pv' r R sg ) - p £ 98 , (15) 

V r d r {pR<j><j>) + (pR<j>4>) V r • V r = ~V • (pv' r R </,</,) - p S 00 (16) 
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v r d r {F r ) + {F r ) V r -v r = 

+P0V(Q)8r 

-(F r }d r v r -(pR rr )d r s 

-(s'd r P') (17) 

-dr (F r v' r ) - 2 (F r v' r ) /r+ (F e v' e + F^ ) /r 
+ (pv' r q/T) , 

and 

pov r d r Q+V r -(p v r )= n „, 
-2(F r )d r s + 2por)(Q.)q/To-V r -(pv' r Q) . us; 

Finally, the spherically symmetric kinetic energy equation is 

v r d r (pK) + {pK)V r -v r = 

+ {p'v' l .)g r +{pR rr )d r v r (19) 
-V,. • (F r ) - V,. • (F P ) + (P'W ■ v') - p Q e . 

3. CHARACTERIZING TURBULENCE OF 2D CORE-COLLAPSE 
SIMULATIONS AND VALIDATING AVERAGED EQUATIONS 

Having derived the mean-field equations and identified the 
important turbulent correlations, we now characterize the 
background and turbulent profiles of 2D simulations. Most 
importantly, we validate that the Reynolds ave rage d equations 
are consistent with the 2D results. Section 13.11 briefly de- 
scribes the 2D simulations, highlighting general qualities that 
are relevant for turbulence analysis such as the location and 
ext ent o f turbulence and neutrino heating and cooling. Then 
in 33.21 we characterize the turbulent correlations in the 2D 
simulations. Finally, in 33.31 we validate the averaged equa- 
tions. 

3.1. 2D Simulations 

The 2D res ults presented here were calculated using 
BETHE-hydro (Mu rphy & Burrowsll2008aj ) and ar e the same 
simulations that were used in lMurphy et al.l (|2009) to develop 
a grav itational wave emiss ion model via turbulent plumes. 
While Mu rphy et al.l (120091) considered a large suite of sim- 
ulations, for clarity, we focus on one simulation that sim- 
ulated the collapse and explosion of a solar met allicity, 15 
Mq progenitor model (IWoosley & Hegerl 120071) and used 
a driving neutrino luminos ity of 3.7 x 10 52 erg s -1 . See 
iMurphv & Burrowsl(l2008al) for more details on the technique 
and iMurphv et al.l (120091) for the setup of this particular 2D 
simulation. 

To demonstrate the evolution of turbulence, most figures 
of this paper highlight three phases after bounce. These three 
stages correspond to modest steady-state convection (404 ms), 
growing convection and SASI (518 ms), and strong convec- 
tion and SASI (632 ms), and the entropy color maps in Fig. 
Q] provide visual context for the shock location, heating and 
cooling, and location and extent of the turbulence. 

Our focus is on the most obvious turbulent region, which 
extends from —80 km to the shock (>180 km). This post- 
shock turbulence is driven by neutrino heating, and in Fig. 
|2] we show neutrino heating (red), cooling (blue), and net 
heating (heating minus cooling, black lines) profiles for ID 
(dashed-lines) and 2D (solid-lines) simulations. These lo- 
cal heatin^_and_£o^ling rates are calculated using eqs. (4-5) 
of IMurphv & Burrows! d2008b) and a neutrino luminosity of 
L Ue = 3.7 x 10 52 erg s" 1 . Below the gain radius, — 100 km, 
cooling dominates heating, but above the gain radius, heating 
dominates cooling. This latter region is called the gain region 
and drives turbulent convection. After matter accretes through 
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Figure 1. Color map of entropy at times representing the three dominant 
phases. The left panel shows the flow during steady-state convection (404 ms 
after bounce), the middle panel illustrates the dynamics and entropy distribu- 
tion during the building convection and SASI stage (5 1 8 ms), and the far right 
panel shows the flow during strong convection and SASI (632 ms). Most of 
the subsequent figures show specific turbulence characteristics at these three 
times. 

the shock, it advects downward through the gain region, pro- 
ducing a negative entropy gradient. In turn, this negative en- 
tropy gradient drives buoyant convection. Though the region 
below the gain radius has a positive entropy gradient and is 
formally stable to conve ction, momentum car ries plumes well 
into the cooling region (IMurphv et al.ll2009l) . This is a well 
known phenomenon in stellar convection and is called over- 
shoot. In neutrino-driven conve ction, the depth of o vershoot 
can be quite large, —20-40 km, ( IMurphv et al.ll2009l) . 

3.2. Turbulent Correlations of 2D Simulations 

Figures[3]and|4]show the radial profiles of the primary cor- 
relations in the averaged equations, eqs. <f4lfT0b. The lowest 
moment correlations are the turbulent enthalpy flux (TqF s ), the 
turbulent kinetic energies (K r & Kg; or Reynolds stresses), and 
the entropy variance (Q), and all three are shown in the top, 
middle, and bottom panels of Fig. [3] Other important higher 
order correlations are the turbulent transport terms, which are 
the transport of entropy flux, (v' r F s ), the turbulent kinetic en- 
ergy flux, (Fk), and the entropy variance flux, (pv' r Q). The 
radial profiles of these higher order correlations are in Fig. [4] 
In general, all turbulent correlations increase over time, and 
the radial profiles indicate that turbulence is dominated by co- 
herent rising and sinking large-scale buoyant plumes. 

The enthalpy (or entropy) flux has two broad physical in- 
terpretations. Most obviously, the entropy flux indicates the 
direction and magnitude of entropy transport due to turbulent 
motions. Naturally, positive and negative entropy fluxes cor- 
respond to upward and downward entropy transport. In addi- 
tion to indicating the direction of entropy transport, the sign 
of the flux also indicates the direction of the buoyancy forces 
driving the turbulence. To understand this second interpreta- 
tion, consider that the entropy flux is defined using the corre- 
lation of the velocity and entropy perturbations (i.e. (s'v 1 )). 
In regions where convection is actively driven, high entropy 
plumes rise buoyantly and low entropy plumes sink. In other 
words, the entropy and velocity perturbations are either both 
positive or negative. Hence, the correlation and entropy flux 
are always positive in regions of actively driven convection. 
At boundaries, where the plumes are decelerated due to a sta- 
ble background, the correlation or entropy flux, is always neg- 
ative. For example, as a low entropy plume penetrates into 
the lower stable layer, the sinking plume becomes immersed 
in a background that has even lower entropy. While the ve- 
locity perturbation remains negative, the entropy perturbation 
changes sign, becoming positive. Consequently, the correla- 
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Figure 2. Neutrino heating (red), cooling (blue), and net heating (black) 
profiles as a function of radius. The times correspond to the stages of post- 
bounce evolution highlighted in Fig. [T] This plot compares the ID (dashed 
lines) and 2D (solid lines) profiles. Above the gain radius (~ 100 km), the net 
heating is positive. In this gain region, a negative entropy gradient formally 
drives buoyant convection. Below the gain region, a positive entropy gradient 
stabilizes against convection. Though only the gain region is formally con- 
vective, overshoot carries convection well into the stable region below (up to 
~40 km below). By definition, the heating profiles are similar. The cool- 
ing profiles show some minor differences, but these lead to minor differences 
in the entropy profiles (see Fig. |6). In Fig. [6] we show that the most no- 
table differences between ID and 2D are produced by the divergence of the 
convective entropy flux. 

tion and entropy flux are negative in the bounding stabilizing 
regions. 

With these interpretations of the entropy (enthalpy) flux, the 
top panel of Fig. [3]shows where convection is actively driven, 
where plumes are decelerated by bounding stable regions, and 
the magnitude of each as a function of time. At the shock, 
th e entropy flux is zero . This is in contrast with the results 
of lYamasaki & Yamadal (120061) . who find an appreciable en- 
thalpy flux at the shock. Whether the enthalpy flux is zero 
at the shock ha s con sequences for the stalled-accretion shock 
solution (See § !4.3l for further discussion). 

Naively, one would expect the gain radius (~400 km) to 
mark the boundary between actively driven convection in the 
heating region and the stabilizing effects in the cooling layers 
below. Though this is roughly correct, careful inspection of 
the enthalpy profiles show that the transition from positive to 
negative entropy flux does not correspond exactly with the 
gain radius. Instead, the gain radius corresponds best with the 
change in slope of the enthalpy profile. Above the gain radius, 



where convection is actively driven, the enthalpy profile has 
a negative gradient, and below the gain radius the gradient is 
positive. 

This profile can be best understood considering a single low 
entropy plume originating at the shock. First of all, as this 
plume accelerates downward, both the entropy and velocity 
perturbations grow in magnitude. This explains the negative 
entropy flux gradient above the gain radius. Below the gain 
radius the background entropy gradient is positive. There- 
fore, the entropy perturbation diminishes as the background 
entropy reduces to the level of the low entropy plume. At 
this point, both the entropy perturbation and entropy flux are 
zero. Consequently, enthalpy flux gradient is positive below 
the gain radius. With a negative gradient above the gain region 
and a positive gradient below it, the enthalpy flux is maximum 
at the gain radius. As the plume's inertia carries it beyond this 
radius, the enthalpy flux becomes negative in the stabilizing 
layers. These general characteristics are observed at all times, 
but with increasing magnitude at later times when convection 
is more vigorous. 

The same simple model can explain the radial ((K r ), green 
lines) and tangential ((Kg), orange lines) kinetic energy pro- 
files in the middle panel of Fig. [3] In fact, like the en- 
thalpy flux, (K r ) peaks at the gain radius. Because the sink- 
ing plumes have higher speeds than the rising plumes and the 
kinetic energy is weighted by the square of the speed, (K r ) 
is dominated by the sinking plumes. Again, the radial pro- 
file of () is consistent with low entropy plumes originating 
at the shock, accelerating downward to a maximum speed at 
the gain radius, and decelerating in the stabilizing region be- 
lo w the gain rad i us. T his is also consistent with the results 
of iMurphv et al.l (120091) . The tangential component, (Kg), on 
the other hand, shows a maximum just 10s of km away from 
the shock. This is where the rising plumes encounter mate- 
rial that has just passed through the shock and both turn their 
trajectories in the ^-direction. 

Interpretation of the entropy variance (Q and the bottom 
panel of Fig. [3]) is less clear. Like the enthalpy flux and kinetic 
energies, Q increases with time. The most naive interpreta- 
tion is that the variation of total heating among the sinking 
and rising plumes increases over time. However, it is not ob- 
vious whether this increase in variance is due to more heating 
or less cooling in either the rising or sinking plumes. ..or both. 
While it seems clear that the profiles of TqF s and K r are domi- 
nated by the sinking plumes, there is circumstantial evidence 
that Q might be dominated by the behavior of rising plumes. 
For one, there is a gradual rise in Q from the lower convective 
boundary to the upper boundary, suggesting growth with the 
rise of buoyant plumes. But the most telling evidence comes 
from the entropy color maps of Fig. [T] In these maps, it is 
obvious that the entropy of the sinking plumes is roughly con- 
stant over time, while the entropy of rising plumes increases 
with time. Hence, the variance of entropy increases with time 
because the maximum entropy of the rising plumes increases. 

All three transport terms in Fig. [4] (v' r F s ), (Fx), and (pv' r Q), 
are negative nearly everywhere. Hence, the flow of core- 
collapse turbulence acts to transport entropy flux, kinetic en- 
ergy, and entropy variance downward. This is typical of buoy- 
ancy driven convection and is observe d in most simulations 
of convection within s tellar interiors (ICattaneo et alJ [l99lT : 
iMeakin & Arnettll2010l) . For the entropy flux, and turbulent 
kinetic energy, this fact further supports the notion that the tur- 
bulent correlations are dominated by sinking plumes. On the 
other hand, at first glance, the negative transport of Q seems 
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Figure 3. Radial profiles for the most fundamental turbulent correlations: 
enthalpy flux, Tq (F s ), (top panel), kinetic energies, K r & Kg, (middle panel), 
and entropy variance, , (bottom panel). The times correspond to the three 
postbounce phases shown in Fig. [T] These plots show that turbulence grows 
with time, and the radial profiles indicate that turbulence is dominated by 
non-local evolution of coherent large-scale buoyant plumes (see text for more 
details). A simple model that considers the evolution of sinking low en- 
tropy plumes explains the featu res of Tq (F s ) and (K r ) and rising high entropy 
plumes explain (Q). See j)3.2l for an explanation. 

to be at odds with our previous conclusion that rising plumes 
dominate the character of Q. However, the moment of the 
entropy variance, (Q), weights only the variance in the en- 
tropy, but the entropy variance flux weights the velocity of the 
plumes as well. In general, the speed of the sinking plumes 
is larger than the speed of rising plumes. Consequently, while 
rising plumes provide the most weight to (Q), sinking plumes 
provide the greatest weight to (pv' r Q). In §|5| we present sim- 
ple models for these transport terms that assume the domi- 
nance of sinking plumes. 

3.3. Comparing Averaged Background Equations with 2D 
Simulations 

Having presented the turbulent correlations of 2D core- 
collapse simulations, we now validate the Reynolds-averaged 
equations. Specifically, we validate the spherically symmetric 
background equations, eqs. (II 1 H 1 3b - For the sake of brevity, 
we do not show a plot of mass conservation but simply re- 
port that M = ATrr 2 {p()V r + (p'v' r )) is indeed satisfied in the 2D 
simulations. 

Figure|5]validates the form of the momentum equation, eq. 
(1121) . including the turbulence terms. In the top panel, we 
plot the velocity profile of ID and 2D simulations as a func- 



Figure 4. Similar to Fig.f5]but for the transport of entropy flux, (v'F s ) (top 
panel), turbulent kinetic energy flux, {Fx), and the entropy variance flux, 
(pv 1 Q) . All transport terms are negative indicating that the transport of each 
correlation is dominated by the large speeds of the sinking plumes. 

tion of radius, and in the bottom panel, we plot the dominant 
force terms in the momentum equation for the 2D simula- 
tions only. Specifically, the bottom panel shows the difference 
in the gravitational and pressure gradient forces, p^g - VPn 
(dashed-line), and the divergence of the Reynolds stress, V • R 
(dotted line). The solid black line shows mV r v r from the 2D 
simulation. This last term is the left-hand side of eq. ( fT2b , 
and in steady state, represents the total force per unit area that 
a Lagrangian parcel of matter experiences. If eq. ( TlZb rep- 
resents the correct derivation of the momentum equation in- 
cluding turbulence terms, then the sum of the right-hand side 
terms (dot-dashed red line) should equal the solid black line. 
Away from the shock, they agree quite well. Interestingly, the 
right-hand side is essentially zero in the heating region where 
convection is actively driven. This implies that the difference 
in the gravitational force and the pressure gradient is nearly 
balanced by the divergence of the Reynolds stress. 

Figure[6]validates the Reynolds-averaged entropy equation, 
eq. ( foi l. The solid black line corresponds to the 2D results, 
and the black dashed line shows the results of ID simulations. 
For comparison, the red curve is the integration of eq. < TT~3T > 
using the ID density, velocity, and heating profiles. Since the 
ID simulations are not able to simulate multi-dimensional ef- 
fects, we omit the turbulent dissipation and the entropy flux 
terms in this integration. The remarkable agreement between 
the results of this integration and the ID simulation bolsters 
our approach in validating eq. ( fT3b . Similarly, the solid green 
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line shows the integration of eq. ( [T3l using the 2D density, 
velocity, and heating curves and no convection terms. This 
curve is similar to the ID results and clearly under predicts 
the entropy in the gain region. On the other hand, including 
the turbulence terms in the integration of eq. dT3b (dot-dashed 
green curve) dramatically improves the comparison. There- 
fore, given the right entropy flux and turbulent dissipation, we 
conclude that eq. ( TT3b accurately determines the background 
flow. Additionally, we conclude that the turbulence models of 
$4]must, at a minimum, produce accurate entropy flux profiles 
to accurately describe the effects of convection in 2D simula- 
tions. 

Figures [3][6] suggest that convection grows monotonically, 
but these figures sparsely sample convection and its effects 
at three times. In Fig. [7J we illustrate that convection in- 
deed grows monotonically from core bounce until explosion. 
To provide some context with shock position, we plot in the 
top panel shock radii as a function of time after bounce. We 
show the ID and average 2D shock radii using dashed and 
solid black lines, respectively. Before the onset of vigorous 
SASI and convection (~550 ms), both stall at ~180 km. Af- 
terward, the 2D average shock radius climbs to ^320 km, at 
which point all measures of the shock unambiguously expand 
in an explosion. To illustrate the asymmetries in the shock we 
also plot the shock radius at the poles (6=0 and 6 = ir). The 
shock radii at the poles oscillate about the average shock po- 
sition until explosion (^700 ms). The middle panel shows the 
maximum of the total, radial, and transverse kinetic energies. 
In general, the turbulent kinetic energy steadily grows until 
explosion. Finally, the bottom panel shows the maximum of 
ID (black) and 2D (green) average entropy profiles. For com- 
parison, we plot the maximum enthalpy flux in orange and the 
maximum entropy resulting from integrating eq. ( TT31 in green. 
As in Fig. |6l the 2D simulation consistently shows higher en- 
tropy at all times, and the results of integrating eq. $13[ are 
consistent with the simulations. 

In conclusion, the results of Figs. (07]) indicate that 
Reynolds decomposition of the hydrodynamics equations, 
eqs. (0]|6]), is consistent with 2D simulations. 

4. TURBULENCE MODELS 



In § |2] we derived the turbulence equations, eqs. (l4lTT0b 
and showed that they suffer from a closure problem. There- 
fore, finding solutions to the background and turbulence equa- 
tions requires a turbulence closure model. In this section, we 
present several turbulence models and identify their strengths 
and weakness es. 

In sections I4.3H4.61 we present four turbulence models. 
These models, with the exception of Model 1, require a model 
for turbulent dissipation and transpor t. Th er efor e, our discus- 
sion of the models begins in sections |4~T1 & 14.21 with a model 
for turbulence dissipation and transport. The next four sec- 
tions, 94. 3114. 51 present models for the primary 2 nd order tur- 
bulent correlations, F s , K, and Q. The fi rst model, Model 1 , is 
a repro duction of a model presented by I Yamasaki & Yamadal 
(120061) and assumes that both (s) and Vso are zero. The re- 
sulting model is simple, but it provides no equation for the 
turbulent kinetic energy. Model 2 is a closure model for 
the Reynolds stress, entropy flux, and entropy variance and 
has been designed and calibrated to simulate isolated buoyant 
plumes. Model 3 is an algebraic model which is akin to MLT 

To model the higher order correlations, Models 2 & 3 use 
expressions involving the local values of 2 nd -order correla- 
tions. Hence, some of the most important terms in a nonlocal 
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Figure 5. Radial velocity (toppanel) and force (bottom panel) profiles for 
the three pha ses s hown in Fig. [T| The forces are the terms in the momentum 
equation, eq. C2J TOP PANEL: The ID velocity profiles for the three phases 
are quite similar. The 2D velocity profiles evolve as the shock radius and 
convective vigor increase. At later times when convection is strongest, the 
velocity gradient tends to zero in the convective region. BOTTOM PA NEL : 
We compare the individual force terms of the momentum equation, eq. (12) . 
The excess force between gravity and the pressure gradient (dashed) is nearly 
balanced by the divergence of the Reynolds stress (dotted). The red dashed 
line shows that the sum of these forces do indeed equal the actual force, m Vv, 
(solid curve). This validates the form of the momentum equation, eq. {12}, 
including the convective term. 

problem are modeled using local approximations. While these 
local models are adequate in some locations, they can be fac- 
tors off in other locations. Rather than relying on these local 
models, we develop in Model 4 (§ 14.61 ) a novel global model 
that uses global conservation laws to constrain the scale of 
convection. 

Later, in § [5] we compare the results of these models to the 
turbulent kinetic energies and entropy flux of the 2D simula- 
tions. 

4. 1 . Turbulent Dissipation: Kolmogorov 's Hypothesis 

Startin g with Kolmogorov 's hypotheses for turbulent dis- 
sipation, lArnett et"aT1 (120091) construct a model for turbulent 
dissipation which they validate with 2D and 3D stellar evo- 
lution calculations. Buoyed by these successes, we construct 
a similar model for turbulent dissipation. One of the primary 
hypotheses of Kolmogorov's theory is that turbulent energy is 
injected at the largest scales and cascades to smaller scales. 
Consequently, the rate of turbulent energy dissipation is gov- 
erned by the largest scales and dimensionally is proportional 
to v' 3 /£, where C is an appropriate length scale, usually the 
largest eddy size. Therefore, our model for turbulent dissipa- 
tion is 

e- — , (20) 
where R = tr(R) is the trace of the Reynolds stress tensor. 
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Figure 6. Radial entropy profiles during the three phases shown in Fig [T] 
In these plots, we compare the results from ID (dashed, black curves) and 
2D (sol id, black curves) simulations to integrations of the entropy equation, 
eq. {13}. The red curve shows the result of this integration using background 
profiles from the ID simulation and no convective terms. The agreement with 
ID simulations validates this technique to verify the equations. The green 
curves show the results of integrating the entropy equation using background 
profiles from the 2D simulations without convective terms (solid green) and 
with convective terms (dashed green) curve. Excluding the convective terms 
under-predicts the entropy profile, and including the convective terms raises 
the entropy and produces accurate entropy profiles. 

In most astrophysical calculations, the length scale is as- 
sumed to be proportional to the p ressure scale heigh t, Hp = 
-{d\nP/dr)-\ On the other hand. lArnett et al.1 (12009b found 
that stellar convection fills the total available space and that 
this also corresponds to the largest eddy size. For very large 
convection zones, they found that C is at most AHp. There- 
fore, we set the length scale as C = max(4f/p,£convX where 
£conv is the size of the region unstable to convection. 

Given Kolmogorov's hypothesis, this model for e is the 
most basic model that one can assume. Later in § 14.61 we 
propose a model for e that satisfies Kolmogorov's hypothesis 
but apparently represents the dissipation of buoyant plumes 
via entrainment. To satisfy global and local constraints for 
convection, we find that e is best modeled by a linear func- 
tion of distance from the shock. In 3D simulations of stel- 
lar convection in which negatively buoyant plumes dominate, 
iMeakin & Arnettl (120101) found a similar result. These results 
suggest that entrainment between rising and sinking plumes 
govern dissipation. 

4.2. Turbulent Flux Models 



Figure 7. A comparison of shock radii and convection diagnostics as a func- 
tion of time after bounce. TOP PANEL: We plot shock radii as a function of 
time after bounce. We show the ID and average 2D shock radii using dashed 
and solid black lines, respectively. Before the onset of vigorous SASI and 
convection (~550 ms), both stall at ~180 km. Afterward, the 2D average 
shock radius begins to climb to ~320 km, at which point all measures of 
shock radii expand outward in an explosion. To illustrate the asymmetries in 
the shock we also plot the shock radius at the poles (9=0 and 9 = it). The 
shock radii at the poles oscillate about the average shock position until ex- 
plosion (~700 ms). MIDDLE PANEL: The middle panel plots the maximum 
of the total, radial, and transverse kinetic energies. BOTTOM PANEL: The 
bottom panel shows the maximum entropy of the I D a nd 2D simulations, 
the maximum entropy resulting from integrating eq. )13t , and the maximum 
enthalpy flux Tq (F s ). Figure[3]samples these convection diagnostics at three 
times and indicates that convection grows monotonically. Theseplots con- 
firm that convection indeed grows with time. Fu rthermore, Fig. [6|validates 
the Reynolds-averaged entropy equation, eq. (13) , at three specific times, and 
the bottom panel validates this equation at all times. 

The gradient diffusion approximatiorQ has found wide 
spread usage for closing the turbulent transport terms (e.g., 
lPopdl2000t Launder & Sandhamll2002l) . For example, Model 
2, which we introduce in 34.41 uses this assumption. The theo- 
retical basis for this closure model is two fold: first, the trans- 
ported quantity behaves like a scalar; second, scale separation 
is a good approximation in the sense that transport is mediated 
by fluctuations on scales small co mpared to the largest scales 
characterizing the turbulent flow (iDaly & Harlowl[l970l) . The 
turbulent transport in a buoyant convection zone, however, is 
generally rec ognized to be mediat ed by larg e scale coherent 
plumes (e.g. ICattaneo et all 1 1 99 It and §13.21 ), thus challeng- 

5 This approximation is also referred to as the down gradient approxima- 
tion, and refers to the closure whereby the flux of a quantity F{ is proportional 
to the gradient of the quantities density Ej, through F{ oc — V£;. 
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ing the theoretical underpinning for a gradient diffusion ap- 
proximation. Furthermore, the shortcomings of the gradient 
diffusion approximation for thermal convection were explic- 
itly illustrated thr ough a series of 3D simul ations of turbulent 
stellar interiors bv lMeakin & Arnettl (12010b . 

Rather than a gradient diffusion approximation, we propose 
flux models that are proportional to R rr x l 2 Ei, where £, is the 
density of the transported quantity. This model is built on 
the advective nature of the transport in which F; oc v'Zs,. We 
propose the following models for turbulent transport of the 
entropy flux, turbulent kinetic energy, and entropy variance 

(v' r F s )n-RH 2 F s , (21) 

(F K )^-p RU\ (22) 

and 

kpv'rQh-^-p^JrQ. (23) 

Except for the entropy variance flux, the comparison with 2D 
simulations ($5]and Fig. |9]l indicates that the constant of pro- 
portionality is ~1. The constant of proportionality for the 
entropy variance flux is found to be ~l/2. 

4.3. Model 1: Steady-State and Zero Entropy Gradient 
Model 

The first turbulence model that we consider is a sim - 
ple model for F s presented by Yamasa ki & Yamadal (120061) . 
lYamasaki & Yam ada (20061) assumed steady state and that the 
entropy gradient in eq. ( fT3] l is zero. From these assumptions, 
they derived a simple differential equation for the enthalpy 
flux. Here, we reproduce this equation, expressing it in terms 
of the entropy flux: 

V,-F, = £. (24) 

A simple integral of this equation with a boundary 
condition leads to an ex pression for the entropy flux. 
lYamasaki & Yamadal ( 12006b assumed zero flux at the lower 
boundary and integrated upward resulting in a nonzero flux at 
the shock. 

A major advantage of this model is that it is simple and 
straightforward to find a solution for F s . Unfortunately, the 
assumptions which lead to simplicity also lead to inaccurate 
turbulent profiles (see § |3}. For one, there is no solution 
for the kinetic energies. Secondly, this model completely ig- 
nores the velocity of the background flow. As we show in 
§ |5] and Fig. Q~2] this is inconsistent with the characteristics 
of neutrino-driven convection in 2D simulations. Thirdly, this 
model produces flawed entropy flux profiles. For example, it 
results in a nonzero entropy flux at the shock, while 2D simu- 
lations show zero entropy flux at the shock (Figure|3]l. In fact, 
the entropy flux is zero at both b oundaries. To compen s ate for 
this nonzero flux at the shock, lYamasaki & Yamadal (12006b 
modified the shock jump conditions. However, the solutions 
to eqs. l4lfT0l are quite sensitive to the form of the boundary 
conditions at the shock. Incorrect boundary conditions will 
lead to erroneous solutions. 

4.4. Model 2: Reynolds Stress and Heat Flux Closure Model 

Model 2 is a Reynolds stress and heat flux closure model 
that has been d esigned and calibrated to model isolated 
buoyant plumes ( lLaunder & S andham 2002). Though core- 
collapse turbulence involves more than just isolated plumes, 



neutrino-driven convection is in fact buoyancy-driven. There- 
fore, it is plausible that the Reynolds stress and heat flux 
model might be an appropriate closure model. 
The model for the Reynolds stress equation, eq. (|7]), is 

V • (pv® R) = p(G + P + n+D - |($ye) , (25) 

where G and P are the buoyant and shear production terms, 
II is the pressure-strain term, D is the diffusion term, and the 
final term is turbulent dissipation. All but the first two terms 
require models, and the expressions for the first two produc- 
tion terms are easily read from eq. 0. 

The first of the modeled terms is the pressure-strain corre- 
lation, and it acts to redistribute energy among the Reynolds 
stress components. For buoyant flows, the pressure-strain cor- 
relation is generally modeled by three contributions 

n = ni+n 2 +n 3 , (26) 

where the first term is proportional turbulent stress, the sec- 
ond term is proportional to the interaction between turbulent 
stress and the mean strain, and the last term is proportional to 
buoyancy. Explicitly, they are 

Ux^-d^Ui-^SijpKj , (27) 

n 2 = -c 2 fp-i%tr(P)) , (28) 

and 

n 3 =-c 3 (G-^, 7 tr(G)) , (29) 

where the parameters c\, c 2 , and c 3 are 3, 0.3, and 0.3 respec- 

ti yely. 

lLaunder & Sandhaml (12002b suggest using the gradient- 
diffusion hypothesis to model the diffusion term. Specifically, 

D = c«V- 0-R-VRJ, (30) 

where the constant c« has been calibrated by experiments and 
simulations to be 0.22. This diffusion term inherently as- 
sumes that the Reynolds stress (or kinetic energy) flux is pro- 
portional to the gradient in the Reynolds stress. This is most 
likely relevant in shear dominated flows. However, the buoy- 
ancy dominated flows of core-collapse convection are charac- 
terized by large scale plumes. As a consequence, the transport 
of kinetic energy flux is directly proportional to these bulk tur- 
bulent motions (i.e. R^' 2 ) rather than the gradient. 

The final term to be mode led is the turbulent dissipation, e. 
lLaunder & Sandhaml (12002b add another differential equation 
that includes more production, diffusion, dissipation terms, 
and constants to be calibrated. To avoid th ese c omplications, 
we simply adopt the turbulence model of § 14.11 

To derive the model turbulent kinetic energy equation, we 
take the trace of the Reynolds stress model equation, eq. d25l ). 

v • V (pK) + (pK) V • v = 

+ <pV)-*-tr((pR>-Vv) (31) 
c s V-(pfR-W)-p e, 

Except for two differences, this model kinetic energy equa- 
tion is very similar to the full kinetic energy equation, eq. ([8]). 
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The first difference is that the work done by turbulent pres- 
sure perturbations, (P'V-v'}, is absent in the model equa- 
tion. Because the pressure-strain correlation, II, is designed 
to only redistribute energy among the components, the trace 
of this term is identically zero. Hence, the model assumes 
that this term is zero. Because low Mach number turbulence 
typically has negligible pressure perturbations, this is a stan- 
dard assumption in turbulence modeling. Our 2D simulations 
confirm this assumption. The second difference is that Fk is 
assumed to be proportional to \7K. This is a consequence of 
the gradient-diffusion approximation, but in § [5] we find that 
Fk is not proportional to the gradient but is best modeled as 
F K ocR 3 / 2 . 

The model for the entropy flux equation, eq. (O, is 

v-VF J +F s (V-v) = 

P0VQg-F s -S7v-R-\7s (32) 
+c,V-(fR-VF,)+n i 

where c s = 0. 15 and the pressure-entropy correlation term (H s ) 
is the analog of the pressure-strain correlation and is modeled 
by three terms 

n s = n sl +n !2 +n !3 . (33) 

In general, these terms act to dissipate F s and represent pure 
turbulence, turbulence and mean strain, and buoyancy inter- 
actions: 

U sl =-c ls ~F s , (34) 
k 

Tl s2 = c 2s F s -Vv, (35) 

and 

n s3 = c 3s p r]Qg , (36) 

where the constants are {cis,C2S,cj,s} = {2.85,0.55,0.55}. 
Once again, the transport term in eq. (f32j) is modeled using 
the gradient-diffusion approximation. 

The final equation models the equation for entropy vari- 
ance, eq. ( [Tol l: 

pv \7Q+QV- (pv) = -2F S ■ Vs c Q V • • Vg) - ^ 

(37) 

where eg = 0.1 1 and r = 0.56. 

Equations d25ll37] i represent a model for the turbu- 
lence equations that is complete and has been exten- 
sively tested and calibrated with experiment and simulations 
(lLaunder & Sand ham 2002j). Unfortunately, there are several 
disadvantages to using this turbulence model. For one, this 
model depends upon a large number of calibrated constants. 
In addition, these equations were designed and calibrated to 
calculate the profile of a single isolated buoyant plume. The 
macroscopic flow of fully developed convection is very dif- 
ferent. Therefore, it is quite possible that the closure relations 
presented in this model are not appropriate for core-collapse 
convection. Furthermore, the dissipation terms make eqs. d25T 
[37l i a stiff numerical problem, requiring careful numerical 
treatment. Consequently, the solutions to these equations are 
extremely sensitive to the uncertain dissipation models. 

4.5. Model 3: An Algebraic Model 

In general , the re are two strategies to finding turbulent so- 
lutions (In § 14.61 we present a third). In the first, the turbulent 
correlations are solutions of differential equations. Models 1 
& 2 are of this type. Within the second strategy, a few key 



assumptions allow the differential equations to be converted 
into a set of algebraic equations, and the turbulence correla- 
tions are solutions to these algebraic equations. An exam- 
ple of a turbulence model which uses a set of algebraic equa- 
tions is the mixing length theory (MLT); it is used extensively 
throu gh out astrophysics, including stellar structure calcula- 
tions dKippenhahn & Wiegerj[T990[) . 

To transform the differential equations into algebraic equa- 
tions, the temporal and spatial derivatives must either vanish 
or be approximated by algebraic expressions. As an exam- 
ple of the first method, one can assume that the temporal and 
advective derivatives are zero, (the L. H. S. of the evolution 
equations eqs. (ITlfTOb). This assumption is valid in most cir- 
cumstances because even though the terms on R.H.S. are large 
they often sum to zero. Therefore, the primary assumptions 
of this approach are steady state and local balancing. In the 
second approach, the spatial derivatives are replaced with ra- 
tios of the variable to be differentiated and a length scale. For 
example, the divergence of the kinetic energy flux is roughly, 
V • Fx ~ Fk/L. In effect, a global boundary value problem is 
reduced to a set of local algebraic expressions. 

In some respects, these algebraic equations still retain non- 
local characteristics. The nonlocality is merely hidden in the 
assumptions. For example, in MLT, local gradients are used to 
calculate the local buoyancy force, but the eddies are assumed 
to remain coherent until a mixing length at which point they 
dissipate their energy. Hence, the finite size of the mixing 
length is an echo of the true nonlocality of turbulence in a lo- 
cal prescription. Where local balancing is important, such as 
the heating region, these approximations can give reasonable 
solutions. However, in regions where nonlocal transport is 
important, such as overshoot regions, these algebraic models 
fail completely. 

To derive an algebraic model, we assume that transport bal- 
ances buoyant driving. Applying this assumption to the exact 
2 nd moment equ ation s, eqs. d8lfT0li. and using the flux models, 
eqs. (I21II23I ). of j34.2l results in approximate kinetic energy, 

-V,-(p ^/ 2 ) ~r)Fsg, (38) 

entropy flux 

-W r -(Rll 2 F s ) ~l poVg Q, (39) 
and entropy variance equations, 

v, • (p„< /2 e) ~4f;Vso. (40) 

The R.H.S. of the approximate entropy flux equation, eq. ( f39b . 
is a sum of two important terms: buoy ancy driving, ppygQ, 
and th e pressure covariance, — (s'W'). lLaunder & Sandhaml 
(120021) suggest that this term is best modeled by -(s'W) ~ 
-(1 /2)porjgQ, and our 2D simulations confirm this model. In 
effect, the sum of these two source terms reduces buoyancy 
driving by a factor of two. 

If we approximate the divergence of a generic flux, Ft, by 
V r -Fj oc—Fi/z, where z = R s -r is the downward distance from 
the shock, then algebraic solutions of the approximate equa- 
tions are proportional to z 1 . To obtain the correct proportion- 
ality constants, we assume solutions of the form R n = az 2 , 
F s = bz 2 , and Q = cz 2 , and substitute these into eqs. ( I38H40I ). 
This results in the following algebraic expressions for the sec- 
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ond moment correlations 

Rrr=\r,g{-d r S )z 2 , (41) 

Fs = ^Po (Tig) 1 ' 2 {-drso? 12 z 2 , (42) 

and 

Q=l(-d r s z) 2 . (43) 

For the entropy gradient in these equations, we assume that 
d r So w Q/(mTo), and to get the radial component of Reynolds 
stress we assume R lr ss R/3 

If instead of assuming a balance between local driving and 
transport, we assume a balance between local driving and dis- 
sipation, then we obtain dimensionally similar results. How- 
ever, they differ in the constants of proportionality and the 
length scale involved. With dissipation, the natural length 
scale is C, which is a fixed scale and is either the full size 
of the convective region (i.e. a constant) or is proportional to 
the local scale height. In either case, using a fixed length scale 
results in algebraic expressions that are inconsistent with our 
2D simulations. Instead, in we show that the convective 
profiles in the driving region of 2D simulations are quite con- 
sistent with z as the length scale. Interestingly, most contem- 
porary formulations of MLT use £ as the appropriate length 
scale, while in the original formulation of the MLT, Prandtl 
used z. In essence, the result that the downward distance, z, is 
the most appropriate length scale suggests that the properties 
of neutrino-driven convection are best modeled by negatively 
buoyant plumes that form at the shock and grow due to buoy- 
ancy. 

4.6. Model 4: A General Global Turbulence Model 

While the previous models are adequate in some locations, 
we show in §[5]that they fail to reproduce the global properties 
of neutrino-driven turbulence. Motivated by this failure, we 
propose an original turbulence model which is derived using 
global considerations. 

Models 2 & 3 employ single-point closure models, which 
assume that on small scales, the higher order correlations can 
be modeled using local 2 nd -order correlations. Consequently, 
some of the most important terms in a nonlocal problem are 
modeled using local approximations. While these local mod- 
els are adequate in some locations, they can be factors off in 
other locations. Given the stiff nature of the governing equa- 
tions, these modest errors can lead to significantly flawed tur- 
bulent profiles. Rather than relying on these local models, we 
use conservation laws and develop an original global turbu- 
lence model. 

Integrating the turbulence equations leads to global con- 
servation laws for turbulence. For example, integrating the 
turbulent kinetic energy equation dictates that global buoyant 
driving equals global dissipation. To satisfy these global con- 
straints, the turbulent correlations relax into the appropriate 
profiles. If the same driving, redistribution, and dissipation 
mechanisms operate in a wide range of conditions, then these 
global constraints also imply self-similar convection profiles. 
For Model 4, we adopt characteristic profiles for the convec- 
tive entropy luminosity and turbulent dissipation. The specific 
profiles are motivated by the evolution of large scale buoy- 
ant plumes and informed by numerical 2 D CCSN simula- 
tions (this paper) and 3D stellar convection dMeakin & Arnettl 



120101) . Then the scales of these profiles are constrained using 
the conservation laws. Given these scaled profiles, we then 
calculate the remaining quantities of interest such as the ki- 
netic energy flux and entropy profile. The latter is particularly 
important in exploring the conditions for successful CCSN 
explosions. 

Our gl obal model builds upon a nd generalizes the ideas put 
forth by iMeakin & Arnettl (1201 (#1 . In the context of stellar 
evolution and no background flow, IMeakin & Arnettl (120101) 
suggest that the primary role of the entropy (enthalpy) flux 
is to redistribute the entropy such that the entropy gradient 
is flat and the entropy generation rate is constant throughout. 
In other words, they found universal profiles for the entropy 
generation rate and entropy profiles. They then integrated the 
entropy and turbulent kinetic energy equations to provide in- 
tegral constraints and solutions for the turbulent scales. 

Similarly, the rate of entropy change in the CCSN context 
is uniform, in fact during the steady-state phase it is zero ev- 
ery where. However, the entropy gradient is nonzero, so the 
IMeakin & Arnettl (120101) model will not suffice. 2D simula- 
tions (see §[5]) indicate that the global constraints and redistri- 
bution drive the convective entropy luminosity, L s = F s AiTr 2 , 
and turbulent dissipation to simple self-similar profiles. 

Informed by 2D simulations of C CSNe (§ El Fig. [131 an d 
3D simulations of stellar convection (IMeakin & Arnettl2010l) . 
we model the convective entropy luminosity with a piecewise 
linear function, 

_ / L s . Q z/zo , forz < zo 

Ls ~ \ L S . Q (Z C - z)/(Z c - zo) , forz > z ' ^ V 

and the turbulent dissipation via 

pe 4irr 2 ~ a(R s -r), (45) 

where z = r-Re is the distance from the bottom of the con- 
vective region, Z c = R S —Re is the total size of the convective 

region, and Re is defined by f^ s Q/Tq dV = 0. In these models, 
the peak of the L s profile is at zo, and goes to zero at Re and the 
shock. The turbulent dissipation is zero at the shock, and in- 
creases linearly downward with a scale of a. Given the many 
differences (i.e dimensionality, accretion, etc.) between 2D 
CCSN and 3D stellar convection simulations, the universality 
of the profiles is remarkable and is a testament to their self- 
similarity. Taken together, this model for L s and e has three 
parameters, L s ,o, zo, and a, which set the scale for convection. 

Next, we use the algebraic expression for the entropy flux 
and global conservation laws to evaluate these scales. Op- 
erating under the assumption that the growth of negatively 
buoyant plumes sets the scale L Si q, we evaluate the algebraic 
expression for the entropy flux, eq. d42b . at the peak: 

L I ,o=47rr 2 J F s , alg (z ). (46) 

We explore two techniques to constrain the position of the 
peak, zo- In the first, we assume that the peak of L„ corre- 
sponds to where entrainment starts to dominate the evolution 
of negatively buoyant plumes and that this distance is propor- 
tional to the pressure scale height (i.e. Z c -zo = ct :o H p ). Com- 
paring to the 2D simulation, we find that a- ~ 1.7, 2.1, and 
2.0 at 404 ms, 518 ms, and 632 ms after bounce. Empirically, 
a Zo ~ 2 with an accuracy of ^10%. While this empirical ap- 

6 These global mod els are similar to the model used for the boundary layer 
in Earth's atmosphere(Plate et al. 1993) 



13 



proach provides a useful diagnostic of zo, it lacks a physical 
derivation. 

In the second technique, we derive a physically motivated 
global constraint for zo- Consider the approximate entropy 
flux and entropy variance equations, eqs. d39l & [40b . They 
represent a set of conservation laws for F s and Q in which the 
source terms are determined by buoyant driving. Combining 
these two equations and integrating over the whole convective 
region, we find that J F s ■ Vs^dr ~ 0. In effect, F s ■ Vso repre- 
sents an important buoyant source term for turbulence, where 
F s • Vso is negative for the growth (decay) of negatively (pos- 
itively) buoyant plumes and vice versa. Because F s = and 
Q = at both boundaries, the integral of this buoyant source 
term must balance out. Therefore, we adjust zo so that F s sat- 
isfies 

F s ■ \7s Q dr = . (47) 



To set a, the turbulent dissipation scale, we first rewrite 
the turbulent kinetic energy equation as an expression for the 
turbulent flux: 

dL k 



dr 



riL s g r - pe4irr 



(48) 



where we have neglected the term due to background advec- 
tion. Integrating eq. d48b over the entire convective region 
and noting that Lk is zero at the boundaries leads to the sec- 
ond conservation law for the balance of total buoyant work 
and dissipation: 



r)L s g r dr-- 



peAirr dr. 



(49) 



Together, eqs. d46ll47l &|49l) constrain the three scales of our 
global model. 

Having set the scales of the turbulent profiles, eqs. d44b & 
(05]), it is now possible to integrate eq. ( 1481 ) to find the turbu- 
lent kinetic flux and, most importantly, evaluate the entropy 
profile including the effects of turbulence. 

Including the time rate of change, the entropy equation is 



d(ps ) . _ Q + pe 
■ +m ■ Vso = 



To 



V • F s , 



(50) 



and integrating this equation over the volume from an arbi- 
trary radius, r up to the shock, R s gives 

|- (J ' ps dV S J+MA So (r)=L s (r)+ £ jrdV+J^ ' jrdV, 

°(51) 

where Aso(r) = so(R s )- so(r) and L s (r) = 4itr 2 F s {r). Assum- 
ing a flat gradient, Asq = 0, and ds^/dt = ((q + e)/7b) m , 
whe re () m is a mass a verage , leads to the integral equation 
that iMeakin & Arnettl (120101) use for stellar convection. For 
the CCSN problem, we assume steady state, dso /dt = and 
A«o 7^0, resulting in 

MAs (r) = L s (r) + J ' ^dV + j " g dV . (52) 

If we momentarily neglect the buoyant term in the kinetic 
energy equation, eq. (1481 ), then the profile for the dissipa- 
tion implies dh^jdr ~ — a(R s — r), which is reminiscent of 
the entrainm ent hypothesi s for the evolution of isol ated buoy- 
ant plumes ( Humeri 1 1973t iRieutord & Zahnlll995l) . In Kol- 
mogorov's hypothesis, dissipation is assumed to be dominated 



by mechanisms at the largest scale. Therefore, eq. (05]) im- 
plies that entrainment of negatively buoyant plumes is the 
mechanism at the largest scales that controls dissipation. 

In summary, we use self-similar profiles of L s and e and 
global conservation laws in Model 4. Assuming that the 
growth of negatively buoyant plumes sets the scale for L s , we 
evaluate the algebraic model, eq. (02]), at the peak (zo) of L s . 
To determine the location of the peak, we set zo such that L s 
and Vso satisfy the global constraint in eq. (07]). Next, the 
dissipation scale, a, is determined by satisfying the balance 
between total buoyant work and dissipation, eq. (l49l . Having 
used global constraints to set the parameters of L s and e, we 
next evaluate the turbulent kinetic luminosity and entropy pro- 
files using eqs. (08]) & (1521 ). Finally, we find R rr by inverting 
our plume model for Fk, eq. 



5. COMPARING TURBULENCE MODELS TO 2D SIMULATIONS 

In this section, we critically compare the results of 2D sim- 
ulations with the turbulence closure models presented in $4] 
Of all the turb ulenc e models that we explore, Model 4, the 
global model (§ 14. 61 , consistently gives the correct scale, pro- 
file, and temporal evolution for the Reynolds stress, R rr , and 
entropy flux, F s . 

Dissipation. — In Fig. [8] we present the time history 
of integrated buoyancy driving and turbulent dissipation. 
The integrated buoyancy work Wt = J (p'v 1 ) ■ g dV should 
balance the integrated turbulent dissipation Qk = J poe dV 
(IChand rasekhar 119611) in steady state. This can be shown 
by integrating the turbulent kinetic energy equation (eq. J8)) 
over the volume of the turbulent region, and assuming that 
the work done by turbulent pressure fluctuations and the 
Reynolds stresses are small, and that the flux of kinetic energy 
is zero out of the turbulent region, all good approximations for 
the scenario under investigation. The volume integrated tur- 
bulent kinetic energy equation then reads 



or 



{p'v')-gdV= / poedV. 



W b = Q K . 



(53) 



(54) 



The overall balance between and Qk presented in Fig. [8] 
shows that the adopted model expression for the dissipation 
(eq. lEOl ) leads to an overall consistency with the evolution 
equation for the turbulent kinetic energy(eq. (8)), and the sim- 
ulation data. 

The global balance between buoyancy driving and turbu- 
lent kinetic energy dissipation tempts one to conclude that all 
buoyancy work is dissipated by turbulent dissipation. While 
this is true in the global sense that the net buoyancy work is 
balanced by turbulent kinetic energy dissipation, it is impor- 
tant to note that the total buoyancy work is an integral over the 
turbulent convection zone of q = (p'v' r ) g r , which is positive in 
the active driving region and negative in the bounding stabi- 
lizing layer. Therefore, some of the work done by buoyancy 
in the active driving region (q > 0) is mitigated by the buoy- 
ancy breaking (q < 0) that takes place in the stabilizing layer. 
It is therefore only the net buoyancy work that is balanced 
by the total turbulent dissipation. During the earliest stages 
buoyancy breaking is significant and is about half the magni- 
tude of dissipation. As convection builds in strength after 250 
ms, the significance of buoyancy breaking steadily diminishes 
until it is about 1/10 of dissipation. 
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Turbulent Fluxes. — In Fig. [9] we present profiles of the 
turbulent flux for three quantities, entropy flux, turbulent ki- 
netic energy, and entropy variance. A comparison to the gra- 
dient diffusion approximation confirms earlier results that it 
is a poor model for transport in thermal convection. Instead, 
we find relatively good agreement with the transport models 
proposed in jj4.2| which are proportional to R rr l ' 2 Ei, where £, 
is the density of the transported quantity. The agreement be- 
tween these models and the data confirm the advective nature 
of the transport (i.e. Fi oc v'Zs,)- 

5.1. Model Comparisons for (Q), {R rr ), and Tq{F s ) 

In this subsection we critically compare each of the turbu- 
lence models introduced in $4] to the simulation data. This 
comparison is summarized by Figures [T0lfT4l The first three 
show the radial profiles of (Q), (R rr ), and Tq (F s ) for the 2D 
simulation data and the turbulence models of § 14. 3114. 51 All 
three figures are divided into three panels with each panel 
showing the 2D simulation profile (solid line) and the turbu- 
lence models at the three times of Fig.Q] The last two figures, 
Figs. [13]&[T4] compare the entropy, entropy flux, and kinetic 
energy flux profiles of the 2D simulation s wi th the profiles 
produced by our global m odel (Model 4, § 14.61 . 

Model 1 co mpa rison ( 34. 3D . — The zero entropy gradient 
(Model 1, § 14.3b model is shown as dot-dashed lines in 
Figs. n~0in~2l Though this model gives the correct order of 
magnitude for (TqF s ), it fails to reproduce the specific radial 
profile and temporal evolution. In fact, there appears to be 
no temporal evolution in this turbulence model, while the 2D 
results clearly grow with time. Furthermore, the zero entropy 
gradient model (Model 1) gives a non-zero entropy flux at 
the shock. This non-zero entropy flux is c learly not a charac- 
teristic of the 2D simulations. See § 14.31 for a discussion of 
how this non-zero entropy flux corrupts the solutions for the 
steady-state accretion sh ock. 

Model 2 comparison ( ? 34.4D . — The Reynolds stress and heat 
flux closure model is represented by a dotted-line in each plot. 
Of all of the turbulence models presented in this paper, it pro- 
duces the poorest reproduction of the turbulence profiles. To 
obtain these profiles we used a shooting method to integrate 
eqs. ( T25ll3~7T i subject to the background flow and boundary 
conditions. At present it is not clear what the specific form for 
the boundary conditions should be, especially at the shock, so 
we used the values of (Q), (R rr ), and (F s ) at a small distance 
from the shock. We then integrated eqs. d25ll37b inward until 
(R rr ) = at the inner boundary. We adjusted the guess for F s 
at the outer boundary so that both F s and R rT are zero at the 
lower boundary. 

For several reasons, we strongly disfavor this model. The 
most obvious reason is the lack of consistency between the 
model and the 2D results. In addition, because these equa- 
tions are stiff, they are quite sensitive to the boundary condi- 
tions and the assumptions for dissipation. We adjusted the dis- 
tance where we sampled the boundary conditions just below 
shock and found the solutions to be extremely sensitive to this 
location. Others have noted similar conver gence problem s 
with boundary conditions to these equations (IWilcoxl 120061) . 
Finally, this model has many parameters that have been cal- 
ibrated for the solution of isolated buoyant plumes, not for 
fully developed convecti on. 

Model 3 comparison ( . ^4.51 ). — In the region where convec- 
tion is being actively driven, the algebraic model, eqs. fill 
l43l , produces reasonably accurate profiles and temporal evo- 



lution. Of the three turbulent moments shown in Figs. [T0lfT2l 
R rr and ToF s matter most in the background equations, eqs. fij- 
|6]), and they show the best correlation with 2D simulations. 
On the other hand, while the algebraic model gives the correct 
scale for the entropy variance, Q, the algebraic model does not 
match exactly the 2D profiles. Fortunately, the entropy vari- 
ance does not directly influence the background equations and 
so in practice this discrepancy can be ignored. Nonetheless, 
this failure should be a clue to what is missing. Furthermore, 
we note that the profiles for R„ and TqF s match the 2D re- 
sults only in the heating region, where convection is actively 
driven. Below the heating region (r < 100km), we set the 
values to zero because it is not clear how to model this region 
with the algebraic model where positive buoyancy de celer ates 
the convective plumes. Finally, as we discuss in § 14.51 the 
success of this model implies that core-collapse turbulence is 
best characterized by low entropy plumes that are initiated at 
the shock, the acceleration of these plumes through the heat- 
ing region, and the deceleration of these plumes at the lower 
boundary by stabilizing gr adie nts. 

Model 4 comparison ( $4.6h — Figures [TTIfT4l show that 
the global model (Model 4) provides the most accurate turbu- 
lent correlations. The Reynolds stress, R rr , and enthalpy flux, 
ToF s , (red solid lines in Figs. [TT1& [T2l derived from Model 
4 have profiles that match the 2D simulation data in scale, 
shape, and temporal evolution. 

Figure [13] compares the terms in the entropy equation, 
eq. 02b . of Model 4 with 2D simulation data, and once again, 
this plot shows that Model 4 accurately reproduces the 2D 
data. The solid blue line represents the change in entropy (in 
units of MAso) due to neutrino heating and cooling alone. 
We find that convection fills the region where this integral 
is greater than zero. In the convective region, the neutrino 
heating and cooling curve accounts for only half of the total 
entropy change at 404 ms and only one third of the entropy 
change at 632 ms. Heating by turbulent dissipation and redis- 
tribution by L s account for the rest. The total entropy change, 
MAso, from Model 4 (red-dashed line) is computed by sum- 
ming the neutrino heating and cooling integral (blue line), 
the modeled turbulent dissipation integral (green-dashed line), 
and the modeled convective entropy luminosity (black-dashed 
line). The modeled entropy difference is only slightly larger 
than the 2D simulation results (red-solid line) and reproduces 
the ge nera l radial profile and temporal evolution. 

In § 14.61 we argue that the global constraints of convection 
and similarity in driving, distribution, and dissipation mech- 
anisms suggest self-similar profiles for L s . Indeed, the cor- 
respondence between our modeled L s (green dot-dashed line) 
and the 2D data (black dashed line) confirms this assumption. 
Moreover, this shape is simply modeled as a piecewise linear, 
pointed hat. The scale of L s is set by the entro py flux we de- 
rive from the algebraic model (Model 3, § 14.51 ) at the position 
of the peak. Since the algebraic model describes the growth 
of negatively buoyant plumes that originate at the shock, the 
scale of L s is in turn set by the growth of these negatively 
buoyant plumes. The position of the peak is determined such 
that the integral constraint, J • Vsodr = 0, eq. iAli . is satis- 
fied. 

In Fig. [14] we compare the kinetic energy flux of Model 4 
(dashed lines) with the results of 2D simulations (solid lines). 
Qualitatively, the modeled fluxes exhibit the correct scales, 
radial profiles and temporal evolution. 

6. TURBULENCE AND CONDITIONS FOR EXPLOSION 
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Figure 8. Buoyancy work integral, W/,, (solid) and total turbulent dissipation, 
Qk (dashed), as a function of time after bounce. At all times, the total buoy- 
ancy work is roughly balanced by the total turbulent dissipation. W/,(q < 0) 
(blue, dot-dashed line) is the buoyancy work for the region where the inte- 
grand of Wt,, q = (p'v' r )g r , is greater than zero, and \Wj,(q < 0) is buoyancy 
breaking where q < 0. At early times, buoyancy breaking contributes sig- 
nificantly to the balance, but is much diminished as the convective motions 
strengthen near explosion. See §|5]for further discussion. 




Figu re 9. C omparison of turbulent transport terms (solid lines) to model 
eqs. )21l23t (dashed lines). This fi gure is similar to Fig. [4] in presentation. 
The transport terms are not propo rtional to the gradient of R rr as the gra - 
dient approximation would suggest (Pope 2000; Launder & Sandham 2002). 

1/2 

Rather, they are proportional to R rr ~Ej, where £, is the quantity in question 
to be transported. This is a consequence of transport by large scale plumes. 

In this section, we suggest that the entropy equation, eq. ©, 
holds the key to understanding the explosion conditions. Fur- 
thermore, we use this equation to argue that of all the con- 
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Figure 10. Radial profile of the entropy variance, (Q) during the three phases 
shown in Fig. [T] In this figure, we compare the results of 2D simulations 
( solid l ine) with the results of the turbulence models. Only Models 2 & 3 
( j|4.4l & l4~5l provide any predictions for (Q), and of these the algebraic model 
(dashed-line, Model 2) gives the appropriate scale and evolution. 



vective terms, the divergence of the convective entropy flux 
most affects the critical luminosity for successful explosions. 
It has been suggested that the critical luminosity condition is 
equivalent to either a ratio of timescales condition dThompsori 
2000tlJankal200UlThompson et alJ20tillMurr>hv & Burrows 
2008b) or to an ante-sonic condition dPejcha & Thompsoii 
201 lb . In either case, it is the entropy equation, eq. (O, that 
leads to this result. For example, if we ignore the turbulence 



terms and integrate the entropy equation over the gain region, 
then we can derive a ratio of advection to heating timescales 



Tadvect 
T heat 



1 



Q_ 

riiT 



dr , 



(55) 



where As = s(r ga ; n ) - s(R s ) is the change in the entropy 
between the shock (R s ) and gain (r ga i n ) radii. Numeri- 
cal results have con firmed that explos i on ensues when this 
ratio exceeds ~1 dBuras et alj 120061: IScheck et al.l 120081: 
iMurphy & Burrowsl2008btlNordhaus et al.l2010l) . In the sim- 
plest interpretation, this result suggests that explosion occurs 
roughly when As(r) exceeds a critical value, As cr i t . Including 
turbulent terms, the change in entropy as function of radius is 



As(r) = - 



R 'Q + Po e, L s {r) 
dr- 



riiT 



M 



(56) 
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Figure 11. Radial profile of the radial part of the R eyn olds stress, (R rr ) dur- 
ing the three phases shown in Fig. [T] As in Fig. \W\ this figure compares 
the results of 2D s imulatio ns (solid lines) with the results of the models pre- 
sented in sections 14.414.61 The 2D results best match the radial profile and 
the temporal evolution of our global model (red line, Model 4). 



where L s (r) = 4irr 2 F s (r) and recall that m and M are negative. 
In our 2D simulations, we find that, of the two turbulent terms 
in eq. d56b , the last term contributes the most entropy change. 
Therefore, if the time condition is a relevant explosion con- 
dition, then the entropy flux is the turbulent correlation that 
most affects the critical luminosity. 

While the timescale co ndition has proven to be a useful 
diagnostic for explosions, iPeicha & Thompson! d201 lb have 
found a more precise explosion condition in that explosions 
occur when c 2 / exceeds ~ 0.2, where c 2 and v 2 are the local 
sound speed and escape velocities squared. They call this new 
condition the ante-sonic condition and note that it varies by 
only a few percent when the neutrino luminosity is changed 
by two orders of magnitude. Over the same range of neutrino 
luminosities, the timescale ratio at explosion varies from 0.7 
to 1.1. In other words, the timescale condition is of order 
1, but it varies by ^1.6. Hence, while the timescale ratio is 
a useful diagnostic relating the most important physical pro- 
cesses, the ante-sonic condition is a more precise (although 
more obscure) condition for explosion. 

Using the integrated form of the entropy equation, eq. 
we show that the timescale diagnostic and ante-sonic condi- 
tion are intimately related. The difference in the sound speed 



Figure 12. Radial profile of the enthalpy flux, Tq (F s ) during the three phases 
shown in Fig. [T] This figure compares the results of 2D simulatio ns (solid 
lines) with the results of the models presented in sections 14.314.61 Of these 
models, only the global model (red line, Model 4) reproduces the scale, pro- 
file, and evolution of the 2D simulation. 

between the shock and an arbitrary radius, r is 



Ac 2 (r) : 



P ci 
po c v 



(57) 



where P and po are evaluated at r, c\ = —(T/p)(dP/dT) p is 
evaluated at constant density, and cy = T(ds/dT)v is the spe- 
cific heat at constant volume. The first term on the R. H. S. 
of eq. (l57l i gives the increase in sound speed due to adiabatic 
compression. The second term represents the change in the 
sound speed due to changes in entropy given by heating and 
cooling, 



and by convection, 



— / -?-dr 

c v mT 



r 2 I 



Cy M 



(58) 



(59) 



Furthermore, since the second term is proportional to the 
change in entropy, it is also proportional to the ratio of 
timescales. Through As, it is apparent that the ante-sonic 
condition, c 2 /v 2 ~ 0.2, is directly related to the timescale di- 
agnostic. 

In a forthcoming paper, we will provide a more thorough 
discussion of these conditions, how they relate to the criti- 
cal luminosity for explosions, and how convection affects all 
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Figure 13. Radial profiles of entropy change in units of MAs. The solid 
lines represent 2D simulation results, and the dashed lines are the results of 
the global model, Model 4. The integrated neutrino heating and cooling (blue 
line) and turbulent dissipation (green dot-dashed line) are evaluated from r to 
the shock radius, R s . The modeled entropy difference (red-dashed line) is a 
sum of the integrated neutrino heating and cooling (blue line), the modeled 
integrated turbulent dissipation (green-dashed line), and the modeled convec- 
tive entropy luminosity, L s (green dot-dashed line). The entropy difference 
derived using the global model (red-dashed line) shows the same scale, ra- 
dial profile, and temporal evolution as the 2D simulation data (red solid line). 
Convection (L s > 0) fills the region where the integrated neutrino heating and 
cooling (blue line) is greater than zero. The 2D simulated L, profile (black- 
dashed line) is self-similar and can be modeled by a piecewise linear, pointed 
hat. 
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Figure 14. Convective kinetic energy flux, Fk, at 404 ms, 518 ms, and 632 
ms after bounce as a function of radius. The results of the global model 
(Model 4, dashed lines) reproduce the scale, general profile, and temporal 
evolution of the 2D data (solid lines). 



three con ditio ns. For now, these analytics reaffirm the suppo- 
sition in § !3.3l that the convective entropy flux most affects the 
explosion conditions. 

6.1. Turbulence and Explosion Condition: 2D Simulation 

Results 

To see if the ante-sonic condition is consistent with our 2D 
simulations, we plot c 2 s /v 2 as a function of radius and at four 
times in Fig. Q3] The first three times correspond to the stages 
shown in Fig. [T] and sample a range of convective strength 
from weakest at the earliest time to strongest at the latest time. 
The final time, 700 ms after bounce, corresponds to the time 
of explosion, which we define as the time when all measures 
of shock radii (see the top panel in Fig. [7]) expand indefinitely. 
For comparison, we show c 2 /v 2 of ID models for these times. 
The ID profiles show very little evolution. However, in 2D 
simulations, the maximum of c 2 / v 2 is strongly correlated with 
the strength of convection and the shock radius. At explosion, 
the peak of c 2 /v 2 is ~0.2, which is consistent wi t h the explo- 
sion condition proposed bv lPeicha & Thompson! d201 ll) . 

It has been argued that convection increases the dwell time 
in t he gain region, which i n turn reduces the critical luminos- 
ity. iPejcha & Thompson! (1201 ll) . on the other hand, propose 
that convection acts to rearrange the flow so that there is less 
cooling, and this reduced cooling is responsible for lower crit- 
ical luminosities. The heating and cooling profiles in Fig. [2] 
and the entropy profiles in Fig. [6] offer a way to investigate 
the merits of each proposal. Unfortunately, interpretation is 
somewhat complicated by the fact that the average 2D cool- 
ing is less than ID cooling for some radii and times but is 
higher for other radii and times. Below ^80 km, 2D cooling 
is always more than ID cooling by ~10%. Above this radius, 
2D cooling is generally a few percent less than ID cooling. 
However, at later times (518 and 632 ms) and above ~420 
km, 2D cooling is a few percent larger than ID again. 

However, Fig. [6] shows that the differences in the average 
cooling between ID and 2D are small and do not greatly affect 
the entropy profile. When the convective terms are ignored, 
the 2D (solid green) and ID (solid red) entropy profiles are 
quite similar. Though there are small differences, the differ- 
ences in average cooling profiles are dwarfed by the effects of 
including the convective terms in the entropy equation (dot- 
dashed green curve). Therefore, it is unlikely that changes in 
the average cooling between ID and 2D lead to the reduction 
in the critical luminosity. Rather, as we argue in $6]it is more 
likely that the divergence of the convective entropy flux is re- 
sponsible for the extra entropy, higher sound speeds, and a 
reduction in the critical luminosity. 

7. CONCLUSION 

Recent simulations of CCSNe indicate that turbulence re- 
duces the critical neutrino luminosity for successful explo- 
sions. This suggests that a theory for successful explosions 
requires a theoretical framework for turbulence and its influ- 
ence on the critical luminosity. In this paper, we develop a 
foundation for this framework, which is represented by the 
following results: 

• We derive the exact steady-state equations for the back- 
ground and turbulent flow. 

• We identify the convective terms that most influence the 
conditions for successful explosions. 

• We have shown that without turbulence, entropy pro- 
files of 2D simulations would be nearly identical to ID 
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Figure 15. Ratio of sound speed and escape velocity squared as a function 
of radius. The four times shown correspond to the three phases in Fig. [T] 
(404, 518, and 632 ms after bounce) and explosion (700 ms after bounce). 
As convection increases in vigor, the maximum in this ratio increases until it 
reaches ~0.2 at the time of expl osion. These results are con sistent with the 
ante-sonic condition proposed by Pejcha & Thompson (201 1). 

and that the convective terms entirely make up the dif- 
ference. 

• This further motivates the need to understand turbu- 
lence in the context of CCSNe. To this end, we cull the 
literature for a broad sample of turbulence models, but 
after a quantitative comparison with 2D simulations, we 
find that none adequately reproduce the global turbulent 
profiles. These single-point models fail because they 
use local closure approximations, even though buoy- 
antly driven turbulence is a global phenomenon. 

• Motivated by the necessity for an alternate approach, 
we propose an original model for turbulence which in- 
corporates global properties of the flow. This global 
model has no free parameters; instead the scale (or pa- 
rameters) of convection are constrained by global con- 
servation laws. Furthermore, this model accurately re- 
produces the turbulence profiles and evolution of 2D 
CCSN simulations. 

Using Reynolds decomposition, we derive steady-state av- 
eraged equations for the background flow and turbulent corre- 
lations, eqs. (fTlfTOl). These equations naturally incorporate ef- 
fects that are important in the CCSN problem such as steady- 
state accretion, neutrino heating and cooling, non-zero en- 
tropy gradients, buoyant driving, turbulent transport, and dis- 
sipation. We validate these equations using 2D CCSN simu- 
lations. For example, we integrate the entropy equation with 
and without the convective terms (see Fig. |6]l. If we neglect 
the turbulence terms, then we recover the ID entropy profile. 
The difference between the ID and 2D entropy profiles is en- 
tirely accounted for by the physics of turbulence. 

Turbulence equations require closure models, but these clo- 
sure models depend upon the macroscopic properties of the 
flow. To derive a closure model that is appropriate for CC- 
SNe, we compare a representative sample of closure models 
in the literature with 2D simulations. Motivated by the fail- 
ure of these models, we have developed an original closure 
model. While the models culled from the literature are single- 
point closure models and use local closure approximations, 
our model is distinguished by using global properties of the 
flow for closure. This global model is further distinguished 
by reproducing the scale, profile, and evolution of turbulence 
in 2D simulations. 



The single-point models use local turbulent correlations to 
derive closure relations for the higher order correlations. Con- 
vection is inherently a global phenomenon, and so while it is 
possible to model the higher-order correlations with local ap- 
proximations in some locations, these models can be factors 
off in other locations. Given the stiff nature of the Reynolds- 
averaged equations, these errors, even if modest, can lead to 
significantly flawed global solutions. Rather than relying on 
these local models, we integrate the turbulence equations and 
derive global constraints based on conservation laws. We pro- 
pose that nonlocal turbulent transport relaxes the turbulent 
profiles to satisfy these global constraints. This relaxation 
combined with the similarity of buoyant driving, entrainment, 
and dissipation leads to self-similar profiles for the most im- 
portant turbulent correlations. In Model 4, we construct a 
global model in which we define these self-similar profiles 
and use global conservation laws to determine their scales. 
Locally, we use the differential form of the conservation equa- 
tions to derive the remaining profiles. 

Our model represents a new approach to turbulence model- 
ing, so we elucidate the assumptions and features that distin- 
guish it from previous models. Single point closure models 
try to employ universal characteristics on the smallest scales 
to close the problem. We are approaching this from the other 
direction. The nonlocal nature of plume dominated convec- 
tion leads us to assume universality on the largest scales and 
a minimum set of global profiles to close the problem. These 
two approaches are complimentary. Assuming universality on 
the smallest scales lends itself to dynamic simulations, while 
the global approach lends itself very well to steady-state prob- 
lems. 

The general strategy that we employ is to establish some 
general characteristic of turbulence and use global conserva- 
tion laws to constrain the scale. For now, we identify the ap- 
parent self-similar profiles as the general characteristic. In 
fact, these self-similar profiles are motivated by the generic 
properties of plume dominated flows and the results of 2D 
CCSN and 3D stellar convection simulations. In the future, 
we hope to identify a more fundamental characteristic and 
physical assumption that leads to these profiles. But until 
then, our global model is the only model that consistently 
gives the correct scale, profile, and temporal evolution for the 
convective kinetic energy flux, Fx, and entropy flux, F s . The 
strongest validation of this model is Fig. [13] in which we re- 
produce the entire entropy profile of 2D simulations. 

In preparation to deriving the reduced critical luminosity, 
we identify the turbulent terms that most influence the con- 
ditions for explosion. Th ree explosion conditions have been 
explored in the literature. iBurrows & Goshvl ( 119931) proposed 
a critical neutrino luminosi ty for successful explosions and 
iMurphv & Burrows! d2008bl) used ID and 2D simulations to 
show that this condition indeed separates steady state accre- 
tion from dynamic explosions. Alternatively, it has been sug- 
gested that explosions occur when the advection timescale 
through t he gain region exceeds the h eating timescale. More 
recently, iPejcha & Thompson] d201 ll) suggest an ante-sonic 
condition in which explosions occur once c 2 s /v 2 e exceeds 0.2. 
In Fig. [15] we show that in 2D simulations, c 2 /v 2 indeed 
reaches 0.2 at explosion. Moreover, using our Reynolds- 
averaged equations, we show that the timescale and the ante- 
sonic conditions are intimately related, and in both conditions, 
convection aides explosion because turbulence raises the en- 
tropy by a term proportional to L s /M. 
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In summary, our global turbulence model contains no 
free parameters, is globally self-consistent, accurately repro- 
duces the mean-field properties of 2D CCSN turbulence, and 
promises to explain the reduction in the critical luminosity. 
Despite these successes, closure approximations generically 
depend upon the properties of the macroscopic flow, mak- 
ing them case-dependent. Hence, it is unclear to what extent 
this turbulence model can accurately describe 3D turbulence, 
especially in the presence of rapid rotation and/or magnetic 
fields. 

Preliminary w ork (iThompson et al.j 12005b 

lYamasaki & Yamadal 12005b iBurrows et al.l 12007 al) hints 
that large rotation rates and magnetic field strengths could 
aid explosion, but it is uncertain how modest values would 
alter turbulence and its effects on explosion. Under the most 
extreme rotation rates and/or magnetic fields, the flow can 
be severely distorted from spherical symmetry (e.g. jets, 
IBurrows et alJl2007ah . In these conditions, it is best to study 
the role of rotation and magnetic fields on turbulence using 
multi-dimensional simulations. On the other hand, for mild 
rotation and magnetic fields, the Reynolds decomposition 
framework employed in this paper can be applied straightfor- 
wardly: mild rotational effects can be i ncluded by retaining 
off-diagonal Reynolds stress terms (e.g.. iGaraud et al.ll2010t) 
and applying Reynolds decomposition to ideal MHD intro- 
duces terms associated with the fluctuations of magnetic 
fields such a s Maxwell stresses and Ohmic heating (e.g., 
IOgilvidr2003h IPessah et all 120061) . However, these analyses 
are beyond the scope of this paper, so for now, we comment 
on the reliability of our global turbulence model for 3D 
CCSN turbulence. 

Even though 2D and 3D turbulence are known to behave 
differently, the global turbulence model reproduces the tur- 
bulent characteristics of 2D CCSN and 3D stellar evolution 
simulations. We suspect, but have not proven, that this is a tes- 
tament to the global nature of the turbulence model. Though 
encouraging, there is no guarantee that the model will work 
so well for 3D CCSN simulations. Therefore, a reliable tur- 
bulence closure model will require comparison with 3D sim- 
ulations. In 2D simulations, steady-state is a valid assump- 
tion. However, differences in the plume structure of 3D turbu- 
lence could lead to more efficient heating, which in turn could 
necessitate including time-dependent terms in the turbulence 
equations. The global nature of turbulence and the similarity 
of driving and dissipation should lead to self-similar profiles 
in both 2D and 3D turbulence. However, the exact profiles 
may differ. Whether any of these differences will affect clo- 
sure approximations is uncertain. Only comparison with 3D 
simulations can clear up this matter. 
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